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Abstract 

A numerical method is developed to solve the time-dependent Dirac equation in cylindrical coordinates 
for 3-D axisymmetric systems. The time evolution is treated by a splitting scheme in coordinate space 
using alternate direction iteration, while the wave function is discretized spatially on a uniform grid. The 
longitudinal coordinate evolution is performed exactly by the method of characteristics while the radial 
coordinates evolution uses Poisson's integral solution, which allows to implement the radial symmetry of the 
wave function. The latter is evaluated on a time staggered mesh by using Hermite polynomial interpolation 
and by performing the integration analytically. The cylindrical coordinate singularity problem at r = is 
circumvented by this method as the integral is well-defined at the origin. The resulting scheme is reminiscent 
of non-standard finite differences. In the last step of the splitting, the remaining equation has a solution 
in terms of a time-ordered exponential, which is approximated to a higher order than the time evolution 
scheme. We study the time evolution of Gaussian wave packets, and we evaluate the eigenstates of hydrogen- 
like systems by using a spectral method. We compare the numerical results to analytical solutions to validate 
the method. In addition, we present three-dimensional simulations of rclativistic laser-matter interactions, 
using the Dirac equation. 

Keywords: Dirac equation, numerical method, cylindrical coordinates, axisymmetric systems, relativistic 
wave packet 



1. Introduction 

The Dirac equation is among the most important equations in theoretical physics and chemistry as it 
gives a quantum relativistic description of fermions such as electrons and quarks. When these particles 
are moving at very high velocity or when they are bound by very strong classical fields, the non-relativistic 
modelling based on the Schrodinger equation fails and theoretical investigations should be based on the Dirac 
equation. The extreme conditions where rclativistic effects arc important can be found in many areas such 
as relativistic heavy ion collisions, heavy ion spectroscopy, cosmology, astrophysics, and more recently, in 
laser-matter interaction (for a review, see [T] and references therein) and condensed matter physics [2J. For 
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this reason, the Dirac equation, coupled to an electromagnetic field, has been studied extensively to evaluate 
many observables such as electron-positron production, molecule spectra, molecular ionization rates, and 
others. However, solving this equation remains a very challenging task because of its intricate matrix 
structure, its unbounded spectrum (the Dirac operator has negative energy states which forbids the use 
of naive minimization numerical methods [3]) and its multiscales (typically, the electromagnetic field is 
macroscopic, with a time scale of £e&m ~ 0.1 — 100 fs = 0.1 — 100 x 10 -15 s, while the electron motion, as 
in the zitterbewegung process, has a time scale of t e i ~ 1 zs = 1CP 21 s, [4]). 

Existing approaches to tackle these important problems can usually be classified in three categories. The 
first one is the analytical method which aims at finding closed-form solutions. Although many important 
problems were treated in this way (HE], it only allows the study of idealized systems. The second approach is 
the semi-classical approximation which can be used to study more complex configurations than the analytical 
method (see [B] for instance). However, it is only valid for a certain range of wave function parameters, 
which may not be realized in the physical system under study. The last one is based on full numerical 
approximations which in principle, can be used to investigate any physical systems. However, even on the 
numerical side, the solution of the Dirac equation is a challenging problem: it requires a lot of computer 
resources [7] and certain numerical schemes are plagued by numerical artifacts such as the fermion doubling 
problem [7J [HI HI HH] • Therefore, special cares have to be taken to resolve these issues when solving the Dirac 
equation numerically for physically relevant systems. 

Among the most successful numerical methods solving the Dirac equation, many are based on a split-step 
scheme where the Dirac Hamiltonian is separated in several operators. This has been used in conjunction 
with spectral schemes in [TTJ HU EH HH US] for the Dirac equation and in [TH] [T7] to solve the coupled 
Maxwell-Dirac system of equations. Very accurate results (with spectral convergence) were obtained with 
these methods. However, the main drawback is that the computation time scales like 0(N log N) (where N 
is the number of spatial points in the discretization) because a Fourier transform has to be computed at every 
timestep. "Real space" methods were exploited by many people using finite element schemes [§1 HE] and 
finite difference schemes (both explicit EH] and implicit [5D] HH H2] ) . However, some of these "real space" 
methods suffer from the fermion doubling problem [HI [23] which induces numerical artifacts and can lead to 
inaccurate solutions. This occurs because the real-time discretization usually modifies the dispersion relation 
[H] such that travelling wave packets acquire a wrong group velocity. Consequently, the phase of a travelling 
wave packet cannot be reproduced accurately by numerical methods suffering from fermion doubling, even 
when the order of convergence is increased |24j . 

Recently, a simple but new numerical method was developed which uses a split-step scheme in "real 
space" and the method of characteristics to evolve the wave function in time, while the space discretization 
is performed with finite volume elements [55] [SB] . In this setting, exact solutions in coordinate space can 
be used in most steps of the splitting (by choosing carefully the time increment St and the element size 
a), resulting in a scheme which is free from the fermion doubling problem and which can be parallelized 
very efficiently |26j . This makes for a very powerful and robust numerical technique which allows to study 
physical systems in Cartesian coordinates, in any number of dimensions. However, for 3-D systems, the 
computational cost is still very important and thus, only short time events can be treated in that case (such 
as heavy ion collisions which last for approximately iRHic ~ 10 -23 — 10 -22 s). For longer events, such as 
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laser-matter interaction (with t pu \ sc ~ 10~ 18 — 10 -13 s), only 2-D calculations are possible and therefore, 
different strategies have to be developed to cope with the high computational requirements. One solution 
is to reduce the 3-D problem to a 2-D problem by using symmetry arguments. In this work, we adopt this 
point of view and study systems which are azimuthally symmetric. For this reason, we extend the split-step 
scheme to solve the Dirac equation in cylindrical coordinates. 

The rationale to consider this coordinate system is twofold. First, many physical systems of interests 
have an azimuthal symmetry and thus, can be treated in cylindrical coordinates. Two examples of this 
are heavy ion collisions at zero impact parameter, and laser-atom interaction in a counterpropagating laser 
configuration. Second, it reduces the mathematical description of a 3-D system to an equation in 2-D 
which of course, reduces the computation time significantly. On the other hand, these coordinates introduce 
new complications in the numerical calculations because the Dirac operator acquires singular terms in the 
coordinate transformation (terms of the form 1/r where r is the radial distance). This complicates the 
numerical evaluation of this operator on the boundary close to r = (dtt r =o,8,z, where Sl r fi. z is the domain 
of the wave function and where r € M + ,6> € [0, 2ir] and z € M). This problem has been studied for other 
equations and many solutions were developed for the Navier-Stokes (and other fluid-like) equations, such 
as the use of pole conditions [57], shifted mesh [25] and series expansion close to the singularity [2S]. A 
treatment of the singularity for the Schrodinger equation in cylindrical coordinates can be found in |30] 
where it is shown that the accuracy of the numerical solution can be improved by writing the differential 
operator in "self-adjoint form". Finally, the Dirac equation with finite difference scheme is treated in [19] 
where a filter is applied at very time-step to get rid of spurious oscillations close to r = 0. In this work, we 
use another approach which consists of a splitting method analogous to the one presented in |25l I26j where 
alternate dimension iteration is performed. The splitting operators are chosen such that all the singular terms 
are included in the radial evolution operator. The resulting equation can then be transformed into a set of 
four 2-D scalar wave equations, expressed in polar coordinates. An integral representation of the solution 
of these equations can then be found: it is the well-known Poisson formula. The latter can be evaluated 
by interpolating the wave function spinor components using Hermite polynomials: using this polynomial 
form allows us to evaluate the integral explicitly. This however entails that a time staggered mesh has to 
be used: the grid points are shifted by a half space step at every time step. This method allows to obtain 
an accurate approximation of the solution at r — while preserving the overall computational performance 
of the numerical method. It also circumvent the singularity problem at r = because Poisson's solution is 
well-defined at that position. Finally, the resulting numerical scheme is very similar to non-standard finite 
difference schemes used to solve other equations [3"T] . 

This article is organized as follows. In Section [2] the numerical method and the discretization of the 
Dirac equation is presented. The splitting scheme is described in details along with boundary conditions at 
r = 0. Section [3] contains several numerical results and benchmark tests. The order of convergence of the 
method is determined numerically by looking at the time evolution of a Gaussian wave packet. The latter 
is also compared to an analytical solution to verify the validity of the method. Then, some more interesting 
physical systems are considered. The first one is a Gaussian wave packet immersed in a counterpropagating 
superintense laser field. It is shown that positive energy states, which can be interpreted as the creation 
of electron-positron pairs, appear in the numerical solution. In the last part of this section, bound state 
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problems are considered where the eigensolutions are determined from an adaptation of the Feit-Fleck method 
[32j to the relativistic case. This allows us to simulate the interaction of a single electron bound in atoms or 
molecules with a strong laser field. We conclude in Section [4] 

Note that in all equations, the light velocity c and fermion mass m are kept explicitly, allowing to adapt 
the method easily to natural or atomic units (a.u.). 



2. Numerical Methods 



The main equation considered in this work is the Dirac equation in cylindrical coordinates which can be 
obtained from the Dirac equation in Cartesian coordinates. The latter is given by [3 3) 



id t ip c (t,x c ) = 



-icd x - eA x (t,x c ) 
-icd z - eA z (t,x c ) 



-icdy - eA y (t,x c ) 



+ fimc 1 + eI 4 V(t, x c ) W) c (t, x c ) 



(1) 



where ip c (t : x c ) E L 2 (R 3 ) <Ei C 4 is the time and coordinate (x c = (x,y,z)) dependent four-spinor, A(t,x c ) 
represents the three space components of the electromagnetic vector potential, V(t,x c ) = Ao(t,x c ) is the 
scalar potential, e is the electric charge (obeying e = — |e| for an electron), I„ is the n by n unit matrix and 
an, j3 are the Dirac matrices. This equation describes physically the relativistic dynamics of a single electron 
subjected to an external electromagnetic field. As usual, the latter is introduced by using the minimal 
coupling prescription^] which allows to preserve the gauge invariance of the equation. 
Throughout this work, the Dirac representation is used where 



The <ji are the usual 2x2 Pauli matrices defined as 



<Ji 
<J l 



(2) 






1 




-i 






, (jy := 




1 





i 



and a z :— 



-1 



(3) 



The Dirac equation in Cartesian coordinates is transformed to cylindrical coordinates by using the fol- 
lowing mapping: 



x 
!) 



r cos 9, 
r sin 8, 



(4) 
(5) 



where the radial distance is r = y^x 2 + y 2 G M + and the polar angle is 9 = arctan(y/a;) £ [0, 2n]. Using the 



The minimal coupling prescription consists of replacing <9 M — > c9 M + ieA^, where is the electromagnetic potential with 
Lorentz index (i. 
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chain rule for derivatives, we gelQ 



+a z 



-icd r — eA r (t, x) 
-icd z - eA x (t,x) 



—ic-dg — eAg(t,x) 



+ /3mc 2 + eI 4 V(i,x) >V{t,x.) 



where x = (r, z, 6) and where the Dirac matrices defined by 



ot x '■= ct x cos a + a y sin ( 
a y := — a x sint 



are now space dependent. It is convenient for the following calculation to write these matrices as 



Ota 

a,. 



a x U(9), 
a y U(9), 



(6) 



(7) 
(8) 



(9) 
(10) 



where U(6) := e eam0ly . When the system has an azimuthal symmetry, the Dirac equation can be simplified 
by using separation of variable. This occurs when has no dependence on the polar angle, i.e. when 
A fl (t,'x.) = Ap(t,r, z). In this case, the Hamiltonian commutes with the angular momentum operator. For 
the remaining of this work, this case will be considered. The ^-dependence can then be factorized by using 
the following ansatz for the four-spinor with cylindrical symmetry |34L 135] : 



^i(i,r,z)e^ 16 
Mt,r,z)e 1 ^ 6 
Mt,r,z)e^ 6 
ip 4 {t,r,z)e^ 26 



U(-6/2)i>(t 1 r 1 z)e ij * f ' 



(11) 



where fj,i t 2 '■= j z T 1/2 and where j z is the angular momentum projection on the z-axis (it can take one of 
the values j z = §,— |, — |, |, |, ...). Inserting this ansatz into Eq. ([6j) and computing the ^-derivative 
we obtain 



iU(-6/2)dti>{t,r, z) 



a x U(6/2) 
+a y U(0/2) 
+a z U(-0/2) 



-icd r — ic- eA T {t, r, z) 



J* 



eAg(t,r,z) 



-icd z - eA z (t,r,z) 



+/3U(-6/2)mc 2 + eU(-6/2)V{t, r, z) \ip(t, r, z) 



(12) 



Finally, using the fact that [U, a z ] = [U,0\ — and a x>y U(6) = U(—9)a XtV , we can completely factorize the 



2 The vector components of the electromagnetic potential transform as A x = A r cos 9 — Ag sin 9 and A y = A r sin 9 + Ag cos 9 
under this coordinate transformation. 
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angle dependence and we get 

id t ip(t,r,z) ■ 



-icd r — ic- eA r (t, r, z) 

2r 



<Xl 



eAg(t, r, z) 



-a. 



-icd z - eA z (t,r, z) 



+ /3mc 2 + eV{t, r, z) >ip(t, r, z) 



(13) 



This equation is the starting point of our analysis and will be solved numerically in subsequent sections. It 
describes physically the wave function for an electron coupled to an electromagnetic field having an azimuthal 
symmetry. 

2.1. Boundary conditions at r = 

In this section, boundary conditions at r — for the wave function are considered (ip(t,r,z)\oQ _ g ). 
When using the cylindrical coordinate system, these boundary conditions are very important to obtain 
physically relevant solutions. In these coordinates, the Dirac equation may have mathematically allowable 
solutions singular at the origin and these have to be discarded as being nonphysicaj^] Moreover, these 
singular solutions can lead to unstable behavior in numerical calculations, in the neighborhood of 9f2 r= o,g,2 



The boundary conditions can be obtained by considering symmetry and regularity constraints of the 
wave function. Following |36j . we make the following assumption: it is assumed that the wave function i\) c 
in Cartesian coordinates is infinitely differentiable (ip c <G C°°) in the whole domain M 3 , and in particular 
for r(x, y) — 0. This assumption, albeit being justified physically, restricts the class of potential that can 
be considered; for instance, the wave function for a Coulomb potential is singular in r = and thus, 
cannot theoretically be considered with our numerical method. Nevertheless, this assumption allows to put 
constraints on the form of the wave function close to r = 0, as shown in the following. 

As demonstrated in [36] , the argument starts with the observation that the coordinate transformation in 
Eqs. Q and §5§ is invariant under the symmetry transformation given by r — >• r' = — r and 9 — > 9' = 9 + n. 
This implies that the wave function in cylindrical coordinates obeys 



This equivalence is shown in Fig. [l] 

Using this result in Eq. ( |TT| ) , we get the following relations on each spinor components: 

Mt,r,z) = (-l)^Mt,~r,z), 

M*>r,z) = (-l) lM2| ^ 2 (t,-r,z), 

Mt,r,z) = (-l)^ 3 (t,-r,z), 

Mt,r,z) = {-l)\^Mt,~r,z). 



(14) 



(15) 
(16) 
(17) 
(18) 



The values /11.2 gives the parity of the wave function with respect to r in the "transformed" polar coordinates. 
Thus, when ^ is even (odd), the components ^1,3 are even (odd) functions of r, while when /j, 2 is even (odd), 



3 This is analogous to the ordinary wave equation for which solutions in the radial direction can be given in terms of Bessel 
and singular Neumann functions. The latter are usually discarded on physical grounds. 



G 



Figure 1: Coordinate transformation symmetry. On the left, the coordinates of an arbitrary point are given 
in usual polar coordinates while on the right, they are given in the "transformed" polar coordinates where 
r£l and 6 G [0,tt]. 



the components ip2,i are even (odd) functions of r. Note here also that if fj-i is even (odd), then /x 2 is odd 
(even); they are not both even or odd at the same time. 

So far, only the symmetry of the coordinate transformation has been used. However, a more explicit 
form for ip can be obtained by utilizing the regularity condition. Using the latter, it is demonstrated in |36j 
that the wave function can be written as 

M*,r,z) = r^h{t,r 2 ,z), (19) 

iP 2 (t,r,z) = r^l/ 2 (t,r 2 ,z), (20) 

M*,r,*) = r^f 3 (t,r 2 ,z), (21) 

M*>r,z) = r^U(t,r 2 ,z). (22) 

where the functions (/»)i=i ... ,4 are regular functions which admit a Taylor expansion in r 2 . These last results 
allow to understand the behavior of the function ip as r — > 0. 

The boundary conditions are given for the following two cases. 

1. When hi is even and fi2 is odd: 

drMtA*) = «> (23) 

MtA*) = 0, (24) 

d r ^{t,Q,z) = 0, (25) 

ip 4 (t,0,z) = 0. (26) 

Moreover, when [Ml 0, we also have tpi(t, 0, z) = ^(t, 0, z) — 0, which implies that both the function 
and its derivative are 0. 
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2. When /ii is odd and [ii is even: 

MtAz) = 0, (27) 

d r Mt,0,z) = 0, (28) 

= 0, (29) 

d r i) 4 (t,0,z) = 0. (30) 

Moreover, when \ii ^ 0, we also have ^(t, 0, z) — ip4(t, 0, z) = 0, which implies that both the function 
and its derivative are 0. 

These boundary conditions guarantee the continuity of the wave equation. They are Robin conditions. The 
parity properties and these boundary conditions will be used in the development of the numerical scheme. 

2.2. Operator Splitting and time evolution 

The main goal of this section is to develop a numerical method that allows to evaluate an approximate 



solution for the Cauchy problem given by a combination of Eq. ( 13 1 with the initial condition 



rl>(t n ,r,z)=,p n (r ) z), (31) 

where the wave function is evaluated at time t n . As previously discussed in 25, 26 for Cartesian coordinates, 
the wave function can then be evolved to a time t n +\ to give ip n+1 (r, z), by using an operator splitting scheme. 
In this work, a similar approach is adopted for cylindrical coordinates and the following splitting is considered 
[2"5l |21>] (note that we omit the r, z in the wave function argument for notational convenience) : 

id t ^\t) = A^(t), ^)(t n ) = tp n , te[t n ,t n+1 ) 
id t ^\t) = B4,^(t), ^ 2 )(t„) = f)(t n+1 ), te[t n ,t n+1 ) 
id t ^(t) = b^){t), ^(t n ) = ^(t n+1 ), te[t n ,t n+1 ) 

and = ^ (3) (<n+l) 

where the upper subscript in parenthesis on the wave function denotes the splitting step number. Note that 
this splitting scheme leads to an error that scales like 0(St 2 ), leading to a first order numerical scheme (for 
more details on the analysis of the method, see [25]). The operators A, B and D will be defined subsequently. 
They are chosen such that an analytical solution can be used to calculate (exactly or approximately) the 
time evolution for every step of the splitting. 

For each step, the initial condition and time domain are shown explicitly. The method consists of solving 
each equation independently with an initial condition given by the solution of the previous step. Note also 
that for every step, the time increment is the same, i.e. 8t n := t n+ i —t n . In the following, as the time step 
is taken constant, we will note St := 5t n 

The operators in the splitting are defined as 

A := —ica x d r — ica x \- ca v — , (33) 

2r r 

B := -ica z d z , (34) 

D := (3mc 2 + el4V(t,r, z) — ea x A r (t,r 7 z) ~ ea y Ag(t,r, z) — ea z A z (t,r, z). (35) 
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This splitting allows to obtain an efficient numerical scheme. For instance, it is shown in the following that 
the second step can be solved exactly using the method of characteristics, using the technique described in 
[25l[26]. The first step, on the other hand, requires a different technique where interpolation and Poisson's 
solution are used. The fulfilment of boundary conditions in r — can be implemented within this framework 
by using an adapted grid and an appropriate interpolation scheme. In the last step, the solution is given by 
a "T-exponential" which is approximated to second order in St. 



2.2.1. Analytical solution for A 

The splitting described in the last section is now used to solve the Dirac equation in cylindrical coordinates 
for systems with azimuthal symmetry. One of the main obstacles to obtain an efficient and accurate numerical 
scheme is the evaluation of the time evolution of the radial operator A. The latter has singular coefficients 
(as 1/r) coming from the cylindrical coordinate transformation. Therefore, the utilization of a naive finite 
difference scheme on the mesh described previously is prohibited as it would require the numerical evaluation 
of A at r = 0, which is obviously undefined. Another approach, suggested in |28j for fluid-like equation, is 
to use a shifted mesh where the grid has no points at r — (the first point is at r — a/2). This method 
however may lead to numerical inaccuracy close to this region in the Dirac equation case. More techniques 
has been developed to solve similar problems in other equations [27J HS1 130] ■ In this work, a new approach 
is developed based on Poisson's integral solution of the scalar wave equation. 



First, the equation obeyed by the wave function is obtained from the definition in Eq. (33) and is given 



by 



idftp^ (t, r, z) = [-ica x d r - ica x ^ + ca y ^] ip^ (t, r, z) 
ik {1) {tn,r,z) = ip n =: g(r,z), 
d t ^(t n ,r,z)=d t iP n =: h(r,z) 



(36) 



Writing the last expression explicitly and decoupling the components allows us to write the following Cauchy 
problem required to evolve the wave function by one timestep, for t € [t n ,t n +i): 



'd?4 1] (t,r,z) 


= c 2 


a r 2 4 


-\Qr~ 


4" 




*), 


%$\t,r,z) 


= c 2 




-\d r - 






z), 


d 2 t ^\t,r,z) 


= c 2 






4 = 




*), 


d^\t,r,z) 


= c 2 






f-'i 
P 


4> ( i\t,r, 


z), 


il>W(tn,r,z) = 


= g(r, 












d t ^(t n ,r, z] 


= h( 


r, z) 











(37) 



The equation obeyed by each wave function component corresponds exactly to the Cauchy problem of the 
2-D wave equation in polar coordinates. An analytical solution to this problem is well-known: it is given 



by Poisson's formula j3jj [38] expressed in polar coordinates (see Appendix B for details on converting 
Poisson's formula from Cartesian to polar coordinates). Using the fact that the angular dependence of the 
wave function can be factorized in the form r, z, 6) = e t/il,2fl ^(t, r, z), the solution shown in Eq. (B.3) 
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yields, for the spinor components i = 1, • • • ,4 



2ncSt J B (r,cSt) 



RdRdO 



yjc 2 5t 2 - [R 2 + r 2 - 2Rrcos(9)] 



x < cos(fi8) 



gi(R, z) + 6thi(R,z) + [R - r cos(6)]d R gi(R, z) 



- sin(^9)—fism(6)gi(R, z) 
R 



(38) 



where fi = [i\ for i = 1,3 and fi = [i 2 for i = 2,4. Also, B(r,c5t) is a disk of radius c5t centered at 
r in the r — (9-plane and corresponds to the integration region (it is depicted in Fig. [4]). As before, 9 is 
the azimuthal angle. The last equation includes the time derivatives hi which are not convenient for the 



hi(r,z) = c 

h 2 {r 1 z) = c 

h 3 (r,z) = c 

h 4 (r,z) = c 



—d r 

-d r 



numerical evaluation. It is possible to discard them by noting that Eq. (36), evaluated at t = t n , gives a 
relation between the different spinor components: 

M2~ 
r - 
Mi" 
r - 

M2" 

r . 
Mi" 

r - 



93(r, z), 
92(r,z), 
9i(r, z). 



Thus, the solution becomes 

ip[ (t n +i,r, z) 



1 



2-KcSt 
x < cos(/^i 



RdRdO 1 
B(r,cSt) ^c 2 5t 2 - [R 2 +r 2 - 2Rrcos(9)} 



51 (i?, z) — cSt 



d R 



M2 

R 



sm(iii6)—^ism(9)gi(R,z) 



W ^n+i, r,z) = 



1 

2nc8t 



g 4 {R, z ) + [R-r cos{d)]d Rgi (R, z) 



1 



RdRdO , 

B(r,cSt) ^c 2 8t 2 - [R 2 +r 2 - 2Rrcos(9)} 



x\ cos(^i 2 0) 



g 2 (R, z) — cSt 



R 



Mi 
R 



g 3 {R, z) + [R -rco$(9)]dRg 2 (R ) z) 



ipf' (t n+ i,r, z) 



s\n{jj, 2 9)—^2 sin(9)g 2 (R, z) \, 



1 



1 



. , RdRd9^= 

2ncSt J B(rtcS t) y/c 2 St 2 - [R 2 + r 2 - 2Rrcos{9)} 



x < cos(mi^) 



g 3 (R,z) - cSt 



Or 



M2 

R 



g 2 (R, z) + [R -rcos{0)}d R g 3 {R, z) 



^4 1) ( < «+l,r, Z) 



sm((j,i6)—(jasm(6)g 3 (Riz) }, 



1 



1 



RdRd9 — 

2itc5t J B (r,cSt) y/c 2 St 2 - [R 2 +r 2 - 2Rrcos(9)} 



x\ cos(^ 2 &) 



94 (R, z) — cSt 



Or 



Ml 
R 



91 {R, z) + [R-r cos{0)}d R g 4 (R, z) 



sin(/z 2 6>) — [i 2 sm(6)g 4 (R,z) }, 



(39) 
(40) 
(41) 
(42) 



(43) 



(44) 



(45) 



(46) 
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which is free from time derivatives. At this point, there is no approximation involved, that is Eqs. (43) to 



( 46 1 are exact solutions of Eq. ( 36 ) . In Section 2.3 the integral in these equations will be evaluated explicitly 
and thus, this solution is at the basis of our discretization method. Finally, is should be noted that if the 
wave function can be Taylor expanded (ip £ C°°nL 2 (IR 3 )<g)(C 4 ), these integrals are well-defined for r £ [0, oo) 
(including the point r = 0) and thus, there is no coordinate singularity problem in this formulation. This 
will be shown in explicit calculation in the following. 

2.2.2. Analytical solution for B 



For the z-coordinate component of the Dirac operator, we have to solve Eq. (32) with B defined in Eq 



(34) on R. To obtain this solution, we starts by diagonalizing the matrix a z using a unitary transformation. 
The resulting system of equations becomes four transport equations which can be solved by using the 
method of characteristics. An inverse transformation brings this solution into the original Dirac matrix 



representation. The details are shown in Appendix A and the final results is given by 



^ 2 \t n+l ,r,z) = U^ + aJ^CWi.n^-c^ + ^-a^^Wi^^ + c^j. (47) 

These solutions represent travelling waves moving in opposite directions at velocity c. A similar technique 
was used in [25l [26] ■ 



2.2.3. Analytical solution for D 

The last equation of the splitting is, for t £ [t n , t n+ i) 

id t ip {3) {t,r,z) = [f3mc 2 + eI 4 V(t,r,z) 

-ea x A r (t,r,z) - ea y A g (t,r,z) - ea z A z (t,r, z)]-0 (3) (t,r, z). (48) 

The solution is simply given by 

■ip ( - i ' ) (t n+1 ,r,z) = fexp^-i[ dT[(3mc 2 - ea x A r (t,r, z) - ea y A e (t,r, z) - ea z A z (t,r, z)] 

ipW(t n+1 ,r,z). (49) 



x exp 



rt„+i 

-ie / oItV(t, r, z) 



where T, the time-ordering operator, has been introduced. The latter orders the argument of the exponential 
function according to their time argument: from the smaller time, on the right to the larger time, on the 
left. This is required because the operator 

h(t) := f3mc 2 — ea x A r (t, r, z) — ea y Ae(t, r, z) — ea z A z (t, r, z) (50) 

does not commute with itself when it is evaluated at different times due to its Dirac structure (the com- 
mutator [h(t), h(t')] ^ 0). In the following, this T-ordering will be omitted (Texp(...) — > exp(...)) as this 
approximation results in an error of 0(5t 3 ) which has the same (or better) accuracy than the method 
considered. It is possible to approximate the time-ordering operator to higher accuracy by using a higher 
order splitting, as discussed in [3"9"ll4"0l HI] . 

It is now possible to discretize the equation spatially by using the analytical solution obtained for all 
operators. 
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2.3. Spatial Discretization 

In this numerical method, the space domain is discretized on a grid forming quadrilateral elements with 
edges of length Sr = 2a and 5z — a (except for elements at r = 0, as described later), where Sr, Sz are 
the distances between grid points in r, z-coordinates and a is the step size. The value of the wave function 
is evaluated on these mesh points and thus, the discretized wave function and electromagnetic field can be 
written as 

Mt,3,k) = P h (j,k)^(t,r,z), (51) 
= ^(t,rf\z fc ), (52) 
Ah, r ,e,z{t,3,k) = Ph(j,k)A rt e,z(t,r,z), (53) 
= A r£ , z {t,rf\z k ), (54) 

where Ph{j,k) is an operator that projects a scalar field on point labelled by (j, k) of the grid h {(j,k) is 
indexing grid points with j £ {1, • • ■ , N r } and k £ {1, • • • , N z } where N r , z are the number of points in r, z 
coordinates, respectively) while iph(t,j,k) and Ah^ r ,e,z(t, j,k) give the value of the discretized wave-function 
and electromagnetic field, respectively, on the grid h. The explicit definition of coordinates r!| and Zfe 
depends clearly on the grid used. In this article, two grids will be considered: hi and h 2 - They are depicted 
in Fig. [2] for the radial direction and the point coordinates are defined by 

a + Sija, (55) 

(56) 

~K (57) 

with z m i n the lower domain boundary in z-coordinates and 8\j the Kronecker delta function. These two grids 
are staggered everywhere except at r — where both have a grid point. This is very convenient because the 
exact value of the wave function or its derivative is known at r — 0, according to the boundary conditions 



Jhi) 
j 


■= U- 


l)2a 


3 


:= (J~ 


l)2a, 




• ^min 


+ {k 



obtained in Eqs. (23 1 to (30). Therefore, the boundary conditions at r — can be implemented exactly on 
both grids. Another remark is that the grid in the z-coordinate is the same for h\ and h^- only the radial 
coordinate grid changes. As explained in the next section, these grids are chosen because the solution of 
the split step in the radial direction requires a time staggered mesh as the mesh changes from hi t 2 to h,2,i, 
respectively, at every time step. We can now discuss the discretization of operators A, B and D on these 
grids. 
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Figure 2: Meshes /ii and h 2 used for the discretization in the radial direction. Circles correspond to grid 
point coordinates r,- . 



2.3.1. Spatial discretization of A 

The discretization of the operator A proceeds by applying the projection operator Ph(j, k) on Eqs. (43) 
to ( |46| and by choosing the lattice spacing such that a — c5t. We obtain 

If „,„,„ 1 

2-na 



n+l, 



,Zk) 



^(t n+ i,rf\z k ) 



,Zk) 



n+l, 



RdRde 



B{rf\a) 



^a? - [B 2 + (rj h) ) 2 - 2Rrf ] cos(6»)] 



x < cos(^ii 



9i(R, Zk) ~ a 



d. 



M2 

R 



g 4 (R, z k ) + [R- rf ] cos(0)]d R gx(R, z k ) 



sm(n 1 e)-^—Hism(0)g 1 (R, z k ) \, 



1 



R 

RdRde - 



(58) 



R? + {rf ) y-2RrV l> cos{0)} 



X \ C0S(^ 2 



g 2 (R, z k ) ~ a 



Mi 
R 



g 3 (R, z k ) + [R- A h) cos(0)}d R g 2 (R, z k ) 



sin(/i 2 0)^M2 s'm(9)g 2 (R, z k ) \, 



R 



(59) 



1 



1 



RdRdO = 

ya 2 - [i? 2 + (rf')> - 2Rrf ] cos(6>)] 



x ^ cos(^ii 



g 3 (R, z k ) - a 



R 



R 



g 2 (R, z k ) + [R- A h) cos(0)}d R g 3 (R, z k ) 



sm(me)-y-nism(0)g 3 (R, z k ) \, 



R 



(60) 



1 



RdRd9 



x I cos([j, 2 t 



g 4 (R, z h ) - a 



Ml 
R 



9i(R, z k ) + [R- rf ] cos(0)]d R g 4 {R, z k ) 



sin(/i 2 ^)^-M2 s'm(0)g 4 (R, z k ) 
H 



(61) 



These equations give a relation between the wave function at t = t n +i and t = t n , assuming that g is defined 



on all points in the interval r £ [rj^ — o J ,f y J a> + a]: this is required to perform the integral on the radial 
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coordinate R over the region B. However, when the equation is discretized, g(R, Zk) is not known everywhere 
on this interval. Rather, we know il>h'(t n ) = gw which is obtained from the preceding time step evolution 
and which is known only on points of the grid h . It is important to notice that a different grid is used for 
this initial condition, denoted by h' ^ h: this choice will be discussed below. Then, we use the following 
strategy to perform the radial and angular integrals: 

• The grid points for b! are chosen as = ± a. Therefore, they are located on the boundaries 
of the integration region B, on the radial axis. 

• An approximation of the wave function g(R,Zk) is obtained by interpolation, such that g(R,Zk) ~ 
gh'{R 1 Zk). Here, g^ is an interpolant of tph'(tn) on the grid h! and thus, can be obtained from 
gh> '■= i>h'(t n )- 



Then, by replacing g — > g in Eqs. (58) to (61 ), it is possible to perform the integrals. 



By making this choice for h' , the two grids at t n+ i and t n are staggered and correspond to hi or h 2 defined 
previously (if h — then h' = /i2,i, respectively). Also, this choice guarantees that only one interpolation 
polynomial is used throughout the integration region B: piecewise polynomial within B is not possible 
because in this case, the integrals cannot be performed analytically. For this reason, the grid points has to 
be positioned on the boundary of the integration region, implying that we have to set cdt = a = Sr/2. In 
the next section, it will be shown that this choice is also very convenient for the longitudinal coordinate 
operator. 

The choice of the interpolation method is also important. First, it has to be a polynomial interpolation 
because in this case, the integral can be evaluated explicitly (this will be shown in the following). In principle, 
any polynomial interpolation could be used, as long as it obeys the following requirements: 



It should obey gt £ C 1 because Eqs. (431 to (46) contain a spatial r derivative 



• It should include information on the derivatives such that boundary conditions at r — are obeyed 
exactly by the polynomial. 

• It should be local such that the implementation is relatively easy to perform and the computational 
performance is preserved. 

• The interpolation error should be bounded by (9(a 3 ) such that the induced error is smaller or at the 
same order as the splitting error. 

In this work, a cubic Hermite interpolation is utilized where the interpolant is given by 

9h'(r,z k ) = g h >(j,k)S 1 (ar))+g h >(3 + l,k)S 2 (t;(r)) 

+(r£2 ~ rf Vfc'(i, k)S 3 (ar)) + (r&'} ~ rf W(j + 1, k)S A (£(r)), (62) 
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Figure 3: Description of the staggered mesh in radial coordinates. The red circle is the integration region 
B. The lines represent the stencil of the scheme. 



(h') 



where £(r) := (hl) 3 (hl) and m/j' is the value of the r derivative of g on grid h! , that is rnh'(j,k) 
Ph'(j, k)d r g(r, z). This interpolation uses third order polynomials given by 

Si(r) = 2r 3 -3r 2 + l, 
S a (r) = -2r 3 + 3r 2 , 
S 3 (r) = r 3 -2r 2 + r, 



S 4 {r) = r 3 -r 2 . 

The derivatives are approximated by a symmetric centered difference as 

9h'(j + l,k) - g h '{j -l,k) 



(63) 
(64) 
(65) 
(66) 



m h > (j,k) 



m h '(j + l,k) 



2(r {h "> -r {h,) ) 
9h'{j + 2,k) -gh'jj, k) 
2(rSS-rf)) 



(67) 
(68) 



It is then a straightforward calculation by using Taylor's approximation to show that g(r, Zk) = gh(f, Zk) + 
0(a 3 ). By using these last two equations, the interpolant can be written as 

9h'(r, z k ) = 9h>(j - 1, fc)Pi(r) + g h >{j, k)P 2 {r) + g h ,(j + 1, fc)P 3 (r) + g h ,(j + 2, fc)P 4 (r), (69) 

where (P n )i<n<4 are third degree polynomials of the form 



p n (r)=£>ro> 



(70) 



(=0 



The coefficients a™(j) can be computed explicitly but are not shown here for simplicity. They depend on 



CO 

J+2- 



The next step is the substitution of Eq. (69) into Eqs. (58) to (61 1, that is we set g(r, Zk) ~ gh>(r, Zk)- 
As a consequence, these equations have integrals of the form 



h{r) 



II rR + 

d6 I dR 



f(6)R l 



a- -/a- -[R- +r- -2Rrcos{9)} 



(71) 
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Figure 4: Integration region B, in blue. The disc is centred on r 
$max, R + and R~~ are the boundary of integration. 



and has a radius of a. The variables 



where f{9) is either cos(/i6*), cos(fiO) cos(8) or sin(/Lt0) sin(0) and where the integer / 6 [0,3]. Also, the 
integration is performed on the disk B(rj : a) centred at rj (the integration region is shown in Fig. |4| and 
the integration on the radial coordinate is performed first. Thus, we have that 

HKVOS ( 1 — ~TT ) , (72) 



2^2 



i?± = r cos(0) ± y/r 2 cos 2 (6>) + a 2 -r 2 . 



(73) 



We now show that these integrals can be performed analytically. First, we make the change of variable given 

by 

R — r cos((9) 



•J r 2 cos 2 (6) + a 2 — r 2 



(74) 



and the integral becomes 

h{r) = 

Now, by using the following result: 

dy 



d0 



i f(0) y^/r 2 cos 2 (0) + a 2 - r 2 + rcos(0) 
dy 



'TT 



_ ri + r-D'i iA^- 
2 1 [ n r(i 



(75) 



(76) 



where T is the Gamma function, it is possible to conclude that all odd powers of the factor yjr 2 cos 2 (9) 



vanish in Eq. (75). As a consequence, after the integration on R, Eqs. (58 1 to (61) only contain terms with 



integer power of cos(#),sin(0),cos(/^6>) and sm(^i9): these can then be integrated easily on angle 9 by using 
well-known equation for trigonometric function integrations. The explicit result of this procedure is not 
shown here for simplicity. It was computed by using the Maple symbolic algebra language and is available 
on request. 



After performing the integrals in Eqs. (58) to (61) by using the interpolant gh>(f,Zk) and the technique 



described previously, the result can then be written in a form reminiscent of a non-standard finite difference 
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scheme (this is depicted in Fig. [3]): 



(h) 



lp2 ) (tn+l,r ( j h) ,Z k ) 



V>3 (*ra+l>''; 



Mj)9i,h(j - 3/2, fc) + A 2 (j)g 1;h (j - 1/2, fc) 
+A 3 {j)gi, h (j + 1/2, fc) + A 4 {j) gi , h (j + 3/2, fc) 
+Bi(j)94,hU ~ 3/2, fc) + B 2 (j)g iih {j - 1/2, fc) 
+B 3 (j)gi,hU + 1/2, fc) + B4(j)94,h(j + 3/2, fc) 
40')52, fc (i - 3/2, fc) + 40> 2 , ft (j - 1/2, fc) 
+40')<72,k(j + 1/2, fc) + A' A (j)g 2ih (j + 3/2, fc) 
+S / 1 (i)53,fe(j - 3/2, fc) + B 2 (j)g 3 ,h(j ~ 1/2, fc) 
+B' 3 (j)93,h(j + 1/2, fc) + B' 4 (j)g 3 M(j + 3/2, fc) 
^i(j')fl3,fc0' ~ 3/2, fc) + A 2 (j)g 3 , h (j - 1/2, fc) 
+A 3 {j)g 3}h (j + 1/2, fc) + A 4 (j)g 3}h (j + 3/2, fc) 
+Bi(j)g2, h (j - 3/2,k)+B 2 (j)g 2ih (j - 1/2, fc) 
+B 3 (j)g 2 , h (j + 1/2, fc) + B 4 (j)g 2 , h (j + 3/2, fc) 
4O')04,fcC?' " 3/2, fc) + A^0> 4 , h (j - 1/2, fc) 
+A' z {])9iAj + I/ 2 , fc ) + A\{j)g ith {j + 3/2, fc) 
+B[(j)g hh (j - 3/2,k)+B' 2 (j)g hh (j - 1/2, fc) 
+B' 3 (j)9i,h(3 + 1/2, fc) + B'S)9i,h{3 + 3/2, fc) 



(77) 



(78) 



(79) 



(80) 



where the finite difference coefficients ^4i,2,3,4(j) depends on the radial position on the mesh. Also, all the 
quantities are written in terms of the grid h such that 



g h (j-3/2,k) = g h ,{j-K_ 2 ,k) 

g h (j -l/2,k) = g h ,{j-K_ x ,k) 

g h (j + 1/2, k) = g h >(j + K u k) 

g h (j + 3/2, fc) = g h ,(j + K 2 ,k) 



(81) 
(82) 
(83) 
(84) 



where the increments are defined as 



if ti = h x < 



K_ 2 = 2 

K-! = 1 

Kx = 
K 2 = 1 



or if h' = h 2 < 



#-2 = 1 
= 

K x = 1 
K 2 = 2 



(85) 



With these results, it is possible to solve numerically the radial part of the splitting involving operator A, 
with an accuracy of 0(a 3 ): the only source of inaccuracy in this case comes from the interpolation error 
because the integrals are computed analytically. The order of convergence of this scheme will be demonstrated 
numerically in Section [3] Note that it could be improved by using higher order Hermite polynomials (quartic 



or quintic). Then, the scheme described in Eqs. (77) to (80) would include more terms 
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The scheme described so far is valid for r > 0. When r = 0, some of the formula are not valid anymore 
and some simplifications can be performed. Moreover, the boundary conditions have to be considered. For 
these reasons, this particular case is treated separately in |Appcndix C[ 

To summarize, the implementation of the numerical scheme consists of the following steps: 

• Evaluate the analytical formula for the coefficients A,A',B,B'. This can be performed by using a 
symbolic mathematical software. 

• Evaluate numerically the coefficients A, A' , B, B' on the grid hi and h 2 . To increase performance, these 
coefficients are stored in a pre-computed table. 



• Evolve the wave function in the splitting scheme by using Eqs. (77) to (80 1. 

2.3.2. Spatial discretization of B 

The spatial discretization of operator B proceeds in the same way as in [33 [5( 
the numerical scheme for this operator by applying the projection operator defined in Eqs. (51 1 and (53) to 



It is possible to obtain 



the analytical solution in Eq. (47). It was shown in [25] that a stable scheme can be obtained by choosing 
a specific relationship between the time and space increments; it should be given by cSt — Ka for any 
K € N*. When this relation is implemented, the numerical solution for this step of the splitting reproduces 
the analytical solution exactly, up to errors related to the projection on the grid and to the splitting, and 
thus, the numerical diffusion is minimized. For K = 1, this is consistent with the previous discussion on the 
radial operator A. 

With this choice of space discretization, the wave function becomes ip^Htn+i, r, z±c5t) —> ip^\t n+ i,j, k+ 
1) and Eq. (47) becomes: 

Wh (*n+l>J)k) = 



4 + o-zM'h (Wi.ii fc - 1) + PU - a z ]tp^ , {t n+1 ,j,k + 1) 



(86) 



Again, we stress that this equation solves Eq. (47) exactly, even after discretization. 



2. 3. 3. Spatial discretization of D 

For the last step involving the operator D, the discretization proceed in the same way as the other 



operator: the projection operator is applied on Eq. (|49|). This yields 



exp 



-if3mc z 5t-ieV£{j, k) 



+iea x A% (j, k) + iea y AZ g (j, k) + iea z A\ (j, k) 



% 2 k ), 



(87) 



U{j, k) exp 



where we defined 



V n (r,z) 



drA rfitZ (T,r, z), 



drV(r, r, z). 



(89) 
(90) 
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The discretization of A™ g z and V n proceeds as in Eq. (51 1 using the projection operator 



The exponential U(J, k) in Eq. (87 1 actually represents a 4 x 4 matrix. It can be evaluated explicitly by 



using a well-known result for the exponential of Dirac matrices, that i^] 



= cos(|F|) + » (/3J? tT F) M\F\) 
\F\ 



(91) 



where F = (Fi, F2,Fs), F, Fi are arbitrary scalar functions and \F\ :— \lF 2 + F • F. Using the last equation, 
the exponential can be given by (here we omit the argument (j, k) in the fields) 



U(j,k) 



c(A)-i*^a(A) 




. A 



Ah ' e+ 2 Ah ' r s{A) 



c(A)-i*^s(A) 



-A ht g+iA ht 



■ A 



s{A) 



-A ht6 +iA htT p 
A 



s{A) 



A h .e+iA h , r 
A 

A 



s(A) 



s(A) c(A) + i^^s(A) 



■ A 



S (A) 



-r^s(A) c{A)+i^s(A) 



(92) 



c(A) := cos(A) , s(A) := sm(A) 



where we defined 
and 

A := yj (mem) 2 + A h>z (j, k) 2 + A h<r (j, k) 2 + A h>e {j, k) 2 . 
This ends the description of the numerical method. 



(93) 



(94) 



2.4- Higher order splitting 

The main numerical error in the scheme presented in the preceding section is due to the operator splitting 
and it can be shown that it is a first order method |25j in time. This can be improved significantly by choosing 
higher order splitting schemes like the second order Strang-like splitting scheme given by (we omit the (r, z) 
in arguments for notational convenience) 





(t) 




^ 1] (t n ) = 


4> n , 


* g [tn,t n+ l) 




(*) 


= B^ 2 )(t), 


^ 2 \t n ) = 




t G [tn,t n+ i) 


idt^ 3. 


(*) 


= DipW(t), 


tf (3) (*») = 




t G [t n , t n+ i) 




(t) 








t G [t n+ i,t n+1 


idtip^ 5 


(t) 


= Aip^(t), 




= * (4) (Wi), 


t G rn+i ! ( n+l 



and V n+1 =V> (5) (Wi) 



(95) 



Here, the time increments are St/2 = t 



n+i 



n+1 



t 



™+5 



for A and B, and finally, St = t 



n+l 



for D. This kind of splitting induces an error of 0(St 3 ). Then, the space discretization proceeds as in the 
lowest order splitting, leading to a second order numerical method |25j . Of course, this can be improved 



'This relation can be derived by using the Taylor expansion of the exponential function and the following Dirac matrices 



properties for any positive integers n: @ 



1, of 



1, /3 2 



2n+l 
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to arbitrary order but this increases the computation time significantly. In this work, we consider only the 
lowest order and second order splitting. 

For an axisymmetric system in 2-D or when the electromagnetic potential has no z-dependence, the 
Dirac equation becomes one dimensional. In this case, the last equation simplifies considerably because the 
number of steps in the splitting is reduced; it is possible to omit every step having the operator B (the same 



is true for the lowest order splitting in Eq. (32)). 



3. Numerical results 

The numerical method exposed in the previous sections is tested in this section by considering some 
physical benchmarks. In the first case, the time evolution of Gaussian wave packets is investigated. This 
is one of the most simple systems and it allows analytical solutions in 2-D. The case of a 3-D wavepacket 
interacting with a laser field is also treated. The other system considered is the calculation of bound states of 
hydrogen-like atoms in 2-D and 3-D. In this case, the eigenenergy is calculated by using a standard method in 
non-relativistic physics: the Feit-Fleck spectral method [31] ■ The latter allows the computation of eigenpairs 
for a given static potential from the solution of the time-dependent equation. The eigenenergies obtained in 
this way are compared to analytical results. 

All the computations are performed in atomic units where the electron mass is given by m = 1.0 and the 
speed of light is c = or 1 where a is the fine structure constant given by a s» 1/137.0359895. Note that in 
these units, we also have the electron charge |e| = 1 and the Plank constant h = 1.0. Also, we use Dirichlet's 
boundary conditions at r = r max , z = z m i n and z max , that is 

dfi r=rmax , 2 , e = 0, (96) 
dtt r , z=Zicim , e = 0, (97) 
dn rtZ=z _ B = 0. (98) 

The value of r max , z m ; n and z max are chosen large enough such that no spurious reflection occurs. Transparent 
boundary conditions are presently under study and will be the main topic of a future publication. 

3.1. Gaussian wave packets and the order of convergence of the radial part 

Gaussian wave packets are certainly among the most elementary systems that can be studied with the 
Dirac equation as they are built from the superposition of plane waves in vacuum (when there is no external 
field). Physically, they represent electrons localized in space, so they have a significant importance in many 
applications. For this reason, they have been studied in many different circumstances [TTl [T2l [TBI H2l H3"I FHI 
FT51 l2l>] . In this study, the time evolution of a simple wave packet will be used to study the convergence and 
to validate the numerical method. Thus, the initial condition considered is given by 

e"^ (99) 



ip(t = 0,r)=N 



r \Ml\ 
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Figure 5: Results for the non-zero components of the Gaussian wave packet with an initial width of A = 0.1, 
evaluated at time t = 0.143 a.u.. The theoretical and calculated curves overlap. 

where Af is a normalization constant and A characterizes the Gaussian width. For \i\ = 0, this wave packet 
represents a spin-up massive electron. The shape of this initial condition is chosen such that conditions 



in Eqs. (19) to (22) are fulfilled at small r while preserving axial symmetry. Then, the numerical scheme 
described in previous sections can be used to evaluate the time evolution of this wave packet. 

When fii = 0, an analytical solution can be calculated (up to a numerical integration) by Fourier 
transform methods, as shown in |Appcndix D| To make a comparison with the numerical results, the 
remaining integrals are calculated in Maple using a numerical scheme based on an adaptive Gauss-Kronrod 
quadrature (with Gauss 30-point and Kronrod 61-point rules), well suited for oscillatory integrands. 

The numerical results are obtained from the scheme described in previous sections, using the first order 
splitting. The width of the wave packet is set to 0.1 and 1.0 a.u. while the domain external boundary is 
''max = 10. a.u., respectively. The time step in both cases is St = 1.113 x 10 -6 a.u.. The non-zero components 
of the wave function are displayed in Figs. [5] and [6] at a time of 0.223 a.u. and 0.557 a.u., respectively. In 
both cases, the calculated results are in very good agreement with the theoretical ones. Also, there are no 
spurious oscillations close to the boundary r — 0, as was observed in [Tj5] for a finite difference discretization. 

Then, the order of convergence of the method in the radial direction is investigated. For the longitudinal 
direction, the numerical scheme is the same as the Cartesian coordinates case, and this was treated in |25| . 
Therefore, the focus here is on the radial operator. 

The methodology used to determine the order of convergence q is based on the evaluation of the numerical 
error in the £ 2 norm defined by 

E#(t) := llVWtW-^WIU, (100) 
where V'cxact is the exact solution. The latter is approximated by the solution calculated on the mesh with 
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Figure 6: Results for the non-zero components of the Gaussian wave packet with an initial width of A = 1.0, 
evaluated at time t = 0.570 a.u.. The theoretical and calculated curves overlap. 



the smallest space step (converged solution). The numerical error is calculated for different mesh sizes and 
angular momenta. The results are shown in Fig. [7] on a logarithmic plot. The slope of the data can be 
easily found by linear interpolation and gives the order of convergence. The numerical values of q are given 
in Table [I] for the considered angular momenta. For all values of j z , the order of convergence in the radial 
direction is close the order of convergence in the longitudinal direction (q 2.0 in the ^-coordinate [25j). 

Table 1: Results for the order of convergence q determined numerically 



jz 


Order of convergence (q) 


1 

2 


1.8999 


3 
2 


1.9435 


5 
2 


1.9767 



3.2. 3-D Gaussian wave packet in a counterpropagating laser field 

One of the important applications of our numerical method concerns the interaction of matter with 
very high intensity lasers. In the last few decades, the laser intensities reached in laboratories have increased 
significantly such that it is now plausible to start observing relativistic and Quantum Electrodynamic (QED) 
effects pp. One of the most important QED observables is the spontaneous production of electron-positron 
pairs from the laser field. Studying this effect from the theoretical point of view requires a solution of the 
time-dependent Dirac equation |221 146| . Therefore, the next example considered is that of a wave packet in 
a time-varying homogeneous electric field. Physically, it corresponds to the interaction of an electron with a 
counterpropagating laser. More specifically, the electric field is given by E(t) = Eof(t) cos(wi) where Eq is 
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Figure 7: Numerical error in the ^ 2 -norm for a wave packet in the radial direction evaluated at time t — 0.0356 
a.u., for different mesh sizes. The lines are linear fits of these errors and the slope gives the order of 
convergence q. 
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the electric field strength, lo is the laser frequency and the envelope function is given by 

-t for t G [0, rnr/u], 

f(t) = { 1 for t e [nir/u, (n + n')ir/u], (101) 

_ {2n+n')^ for te ^ n + n ') 7r / Wj ( 2n + n')n/uj} 7 

where n, n' £ N counts the number of half-cycles for linear ramping and for a constant enveloppe, respectively. 
Note that there is no space dependence in E(t), which corresponds to the field of the laser at the anti-node 
of the standing wave. We work in a gauge where A n = and thus, the electric vector potential is given by 
A x (t) = / dt'E(t'). The initial state represents a positive energy wave packet at rest and is assumed to be 



given by Eq. (99) with /ii = (we consider j z — 1/2). 

The simulation is shown in Fig. [8] for cu = 100 a.u. and E Q = 3.65 x 10 6 a.u., which is slightly higher 
than Schwinger's critical field (2.6 x 10 6 a.u.) at which static electron-positron pair production starts to be 
important. The width of the initial wave packet is set to 1.0 a.u. while the domain external boundaries are 
at ?*max = 10. a.u., 2 m i n = —20. a.u. and z max = 20. a.u.. The time step is St = 4.86 x 10 -5 a.u.. 

The numerical results show that secondary peaks are formed at each cycles. The group velocities of 
these peaks and that of the primary peaks are in opposite direction and thus, they are carrying different 
electric charge. This can then be interpreted as the production of electron-positron pairs [461 113) . This 
is because we started with an initial state containing only positive energy states. By interacting with the 
intense electromagnetic field, negative energy states are created. It should be noted however that a full 
calculation of the pair production rate requires a summation over all the positive and negative energy states 
[4"6"j : this is not performed here. 
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Figure 8: Time evolution of the density IV'I 2 for a wave packet in a counterpropagating laser field at (a) 
t = 0.0 a.u., (b) t = 0.044 a.u., (c) F = 0.088 a.u. and (d) t = 0.131 a.u.. 
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3.3. Bound states 

It is very important in many applications to be able to evaluate bound states of the Dirac operator 
because many observables can be related to these entities. Also, the static bound states can serve as initial 
conditions for the dynamical evolution of the system under study. For these reasons, many numerical 
methods were developed to compute the Dirac eigenenergies and eigenpairs for a given static potential 
[551 |4"TI |4"51 liUl I5UI 15T1 I5"2"] . It is possible to compute time-independent wave functions from a time-dependent 
numerical scheme by using the well-known Feit-Fleck method [35], as explained in the following. 

3.3.1. Spectral method 

In this section, the spectral method used to compute the eigenenergies is presented. Most of this section 
is based on [52] where this method was used to evaluate the eigenstates of the Schrodinger equation and 
where more details can be found. This method was also used for the Dirac equation in JT5J [T5J [S5J D2] to 
calculate the eigenfunction for hydrogen-like atoms. 

The main ingredient of this numerical scheme is the auto-correlation function defined by 

C(t) := (*(0)|*(t)>, (102) 
cFx.&(Q,x)V(t,x). (103) 



where ^(0,x) is an arbitrary trial function. When the potential is static (independent of time), the wave 
function can be expressed as a superposition of eigenstates. It can then be easily demonstrated that the 
Fourier transform of the auto-correlation function C(E) is sharply peaked at the eigenenergy values; C(E) 
becomes a sum over Dirac delta functions positioned at the exact bound state energies. This however 
implies that the wave function should be evolved to an infinite time to evaluate the Fourier transform, which 
of course, is impossible in a numerical calculation. 

The same procedure can however be performed when only a finite time is available. In this case, the 
Fourier transform is calculated as 

/oo 
dtw(t)C(t)e iEt (104) 
-oo 

where w(t) is a window function. The latter allows to determine the functional form of the lineshape, i.e. 
the Fourier transform of the window function gives the equation of the eigenenergy peaks. It is convenient 
to choose the Hann window function given by 

f fortG [0,71 

w(t) = \ , (105) 

[0 forte (-oo,0) U(T,oo) 

where T is the final time of the calculation. This choice of window function allows an accurate determination 
of resonance position because the resulting lineshape has no side lobes which could be confounded with other 
resonances. 

As seen previously, the eigenenergies can be determined by looking at the power spectrum of the auto- 
correlation function. Once the energy E of the bound state is determined, the corresponding eigenstate 
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can be calculated by using 

* E (x)=/ dW{t,x)w{ty E \ (106) 
Jo 

where ^(OjX) is an arbitrary trial function. 

The numerical procedure can be summarized as follow: 

1. Determine the eigenenergy 

• Choose a trial function. 

• Evolve the trial function in time by computing C(t) at each time step. 

• Compute the Fast Fourier Transform (FFT) of C(t)w(t) to obtain the power spectrum. 

2. Determine the eigenfunction. 

• Choose a trial function 

• Set the value of the energy to the one determined in the previous step. 



• Evolve the trial function in time by calculating Eq. ( 106 ) at each time step. 

This numerical method requires at least two calculations of the time-dependent Dirac equation to obtain the 
eigenstate. This is clearly not as efficient as other methods such as variational schemes [3l|48l[50]. However, 
it is very easy to implement and it allows to evaluate the eigenstate of the grid used in the split-step method; 
in contradistinction with basis set expansion methods which can not be easily adapted to our numerical 
method. 

3.3.2. 1-D exponential potential 

By combining the Feit-Fleck method with the scheme for solving the time-dependent Dirac equation, it 
is possible to evaluate the spectrum and eigenfunction of any bounding potentials consistent with the chosen 
grid. The first test of the spectral method is for a simple exponential potential in 1-D. We consider a static 
scalar potential given by 

V{z) = -^e-^ (107) 
la 

where g and a are parameters characterizing its strength and width, respectively. The Dirac equation can 
be solved exactly in this case and the eigenenergies E are solutions of (in units where c = 1)|54) 

E — iq\ 2 \Fi(qa + iEa, 1 + 2qa; ig)\F\ (1 + qa — iEa, 1 + 2qa\ — ig) 



\Fi(qa — iEa, 1 + 2qa; — ig)iFi(l + qa + iEa, 1 + 2qa; ig) ^ ^ 

where q = \J m 2 — E 2 and i-Fi is the confluent hypergeometric function. This transcendental equation can be 
solved numerically and thus, can be compared to the numerical results obtained from the spectral method. 
For calculation using the Feit-Fleck method, we start with a trial function given by a gaussian wave 



packet representing a spin up static electron (see Eq. (D.l), but set x = z and y = 0). We consider a 
1-D domain, the z-axis, and make 2 x 10 6 time iterations to get the required precision. The power spectra 
obtained for two cases (g = 1.0; a = 2.0 and g = 2.0; a = 1.0) are shown in Fig. [9] We find ground 
state energies of -E gr0 und = 16875 and Aground = 11218, for the two cases, respectively, which is close to the 
analytical solution, where we get E = 16866.06 a.u. and E = 11201.44 (the relative difference is less than 
0.15%) 
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Figure 9: Power spectrum for the exponential potential for two cases: g — 1.0; a = 2.0 and g = 2.0; a = 1.0. 
The peak is at the position of the ground state energy. The analytical solution yields ground states of 
E = 16866.06 a.u. and E = 11201.44 a.u., respectively, so the peaks are at the right position. 



3.3.3. 3-D hydrogen-like atom in a laser field 

This last example concerns the interaction of a laser field with the ground state electron of a hydrogen-like 
atom (Is orbital). The ion is modelled by the 3-D Coulomb potential. The latter can be solved analytically, 
when no laser field is present, and represents physically the electric potential of a point charge. Numerically 
however, it is problematic because the wave function and potential have a singularity at the charge position, 
in r = 0. Thus, it does not obey some of the assumptions required by the numerical method presented in this 
article. For this reason, rather than using the Coulomb potential directly, a regularized potential is used. It 
is given by 

{tr(w-Z) for r<R 
V{r) = } 2R \ R J (109) 

[ f forr>i? 

where R here is the radius of the nucleus and Z is the atomic number. For r < i?, this potential corresponds 
to the electric field of a charged sphere having a radius R with constant distribution of charge. In this work, 
this is chosen to insure that ip £ C°° such that relevant numerical results can be obtained. In the following 
calculation, the radius of the nucleus is set to R = 0.01 a.u. while the atomic number is Z = 10. 



The initial state is prepared by using the spectral method presented in Sec. 3.3.1 The trial state is a 
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Figure 10: Results for the power spectrum of trial wave function and the bound state constructed from the 
Feit-Fleck method. The Feit-Fleck method filters the wave function such that only the bound state remains. 



wave packet centred on the atom, it is given by 



ip(t = Q,r,z) = N 



r \p-i\ 







e 4 ^ 2 . 



(110) 



In the first step of the calculation, where the spectrum is evaluated, this initial state is evolved up to 
tj = 47.5 a.u., making for an energy resolution of 8E ss 0.13 a.u. in the spectral method. It was then 
determined, from the power spectrum of the trial wave function shown in Fig. |10| that the ground state 
energy is -E gr0 und ~ 18710.3 a.u.. This value is close to the analytical Coulomb ground state energy given by 



18729.9, with a relative difference of 5 re \ « 0.1%. The ground state is then constructed from the 



same trial function by using Eq. (1061. The spectrum of this bound state is also evaluated as a consistency 



check and is depicted in Fig. [10] A comparison of the trial state and bound state spectra demonstrates 
clearly that the spectral method filters out the unwanted frequencies while keeping only the ground state 
component. The resulting ground state is shown in Fig. |11| 

The wave function representing the ground state of an hydrogen-like atom is then evolved in the field of 



a counterpropagating laser pulse, for which the vector potential is given in Eq. (101). The laser parameters 
are the same as in Section 3.2 except for the maximum electric field which is set to E = 3.65 x 10 5 a.u.. 



The evolution of the wave function is shown in Fig. 12 It can be seen in these pictures that the electron 
moves from left to right, driven by the intense laser field. 
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Figure 11: Ground state of a hydrogen-like atom with Z = 10. 



4. Conclusion 

In this work, a numerical method based on a splitting scheme was developed to solve the time-dependent 
Dirac equation in cylindrical coordinates. The main new idea was the utilization of the Poisson integral 
solution and Hermite interpolation to solve the radial part of the splitting. This allowed us to circumvent 
the coordinate singularity problem in r = 0, while performing all the calculation in coordinate space (no 
Fourier or Hankel transforms are required): the resulting algorithm is very similar to a non-standard finite 
difference scheme. 

This numerical scheme is an extension of the one presented in [251 126] and thus, shares some of its main 
properties: 

1. There is no fermion-doubling problem in the discretization process. Indeed, each step of the splitting 
is solved in a way that does not modify drastically the continuum dispersion relation, apart from small 
errors related to the splitting and the radial spatial discretization. A careful analysis of this issue 
would require a Von Neumann analysis, as in |26j for Cartesian coordinates, where the continuous 
and discretized Dirac equation are Fourier transformed^ to obtain their respective dispersion relations 
(similar calculations are found in [23 HI])- This will be the subject of future investigation. 

2. The numerical scheme can be parallelized very efficiently using a domain decomposition method be- 
cause the algorithm obtained is local in space. Thus, the domain can be separated into subdomains 
which can be solved independently on different processors: only the information at the subdomain 
boundaries is exchanged between neighbour processes. This leads to a quasi-linear speedup as the 
number of processors is increased |26) . This was not considered explicitly in this article, but was 
verified numerically in some calculation. 

The method was tested by looking at the dynamics of wave packets in 2-D. It was shown that the 
numerical results were in agreement with the analytical result. Moreover, it was tested that the order of 
convergence of the method in the radial direction is similar to the Cartesian case. Then, a more interesting 



5 For the radial part, we should use the Hankel transform. 
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Figure 12: Time evolution of the density \ip\ 2 for an electron in the potential of an hydrogen-like atom (with 
Z = 10) in a counterpropagating laser field. 
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system was studied which required much more computational resources and where the efficient parallclization 
of the numerical scheme was required. This allowed us to simulate a 3-D wave packet in a counterpropagating 
laser field and show the appearance of both negative and positive energy states (interpreted as the creation 
of antimatter) during the time evolution. 

The numerical method was also combined with the Feit-Fleck spectral scheme and allowed us to evaluate 
the eigenenergies and eigenfunctions of a 3-D Coulomb-like potential representing a nucleus. Although 
the performance of the Dirac solver developed in this paper is very good, the calculation of the eigenpairs 
required a massive amount of computation time. This is mainly due to the slow convergence of the Feit-Fleck 
method for which the energy resolution scales with the simulation time: a higher resolution requires more 
calculation times. For this reason, it would be very challenging to make 3-D Quantum Electrodynamics 
(QED) calculations with the combination of these two methods as the latter necessitates a summation over 
intermediate states, which involves the computation of all eigenfunctions. However, it may be possible to 
use a different method to evaluate the initial state which converges faster than the Feit-Fleck method. This 
is presently under investigation. 
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Appendix A. Solution of the Dirac Equation in z-coordinates 

This Appendix describes precisely how the solution of the Dirac equation in z-coordinate is obtained. 
We are interested first in solving the following equation: 

id t ip(t, r, z) = -ica z d z ip(t, r, x) (A.l) 

with an initial condition given by 

i>(t ,r,z)=g(r,z). (A.2) 

In the first step, the matrix a z is diagonalized to decouple the spinor components. This is performed by a 
similarity transformation as 

P\a z P z = A z (A.3) 

where P z is a unitary matrix and A z = Diag[l, 1, —1, —1] = is a diagonal matrix. Starting with a z in the 
Dirac representation, the explicit expression of the transformation matrix is 

P z = ^(P + a z ). (A.4) 

The resulting Dirac equation is then 

id t ip{t, r, z) = -ic/3d z ip(t, r, z) (A. 5) 
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where we defined ip := P^tp. It is convenient here, for notational purposes, to split the four-spinor into two 
bi-spinors as i\> = (</5, x) T to get 

idt<p(t,r,z) = —icd z ip(t,r, z) (A. 6) 

idtx(t,r,z) = icd z x(t,r,z) (A. 7) 

Therefore, the Dirac equation clearly becomes a set of four uncoupled first-order differential equations in 
this representation. Their solution is well-known and can be obtained from the method of characteristics. 
The solution is given by 

<fi,2(t,r,z) = g lt2 (r,z - cSt), (A. 8) 

Xi,2(t,r,z) = g 3A (r,z + c5t), (A. 9) 

along with the conditions z ± zq = cSt. The latter is the characteristics along which the partial differential 
equation becomes an ordinary differential equation. Note that the initial conditions are related to the original 
representation as g — P\g. 

In the last step, the solution is transformed back to the original representation. The final result, after 
some basic manipulations, is given by 

ip{t,r,z) = -{[I4 + a z ]g(r,z-c5t) + \l4-a z ]g(r,z + c5t)}, (A.10) 

This equation is used in the numerical method to evolve the wave function in time in alternate direction 
iteration. 

Appendix B. Solution of the 2D Wave equation 

The Cauchy problem for the 2-D wave equation in Cartesian coordinates is given by 

d t u(x,t) = c 2 V 2 u(x,i) 

u(0,x):= Uo (x), (B.l) 
d t u{0,t) := u (x). 

for a general regular u. It is well-known that the solution to this 2-D problem is given by the Poisson formula 
[371 [35]: 

«(x,t) = ~ I d y 1 =My) + fa (y) + Vuo(y)-(y-x)] (b.2) 

27rci J B ( x ,ct) yjc 2 t 2 - |y-x| 2 

where -B(x, ct) is the ball of radius ct about the position x. Writing this equation in cylindrical coordinates 
and assuming that the solution describes an axisymmetric system such that the angular dependence can be 
factorized as u(r,9,t) — e lf * e u(r 7 t), we obtain 



If 1 

u(r,t) = / RdRdO 

2nCt J B (r,ct) 



/s(r,ct) ^cH 2 - [R 2 + r 2 - 2Rr cos(6>)] 

x \ cos(/if) u (R) + tv a (R) + [R - r cos(e)]d R u (R) 



sm(fj,6)^-iJ,sm(e)uo(R)\. (B.3) 
R 
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This gives an integral representation for the solution in cylindrical coordinates of the wave equation given 
by: 

,2 ~ 



d t u(r,t) = (? [d 2 r + -d r 



u(r, t). 



(B.4) 



This equation is the same as those found in Eq. (37) and therefore, their solutions are given by Eq. (B.3|. 



Appendix C. Spatial discretization of A when r — 



When r = 0, Eqs. (58) to (61) become 



tp[^ (t n+1 ,0,z k ) 



i r i 

/ RdR -. , 

2na J yja 2 - R 2 Jo 



gi{R, z k ) - a 
RdR 
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R 



d6 cos(/ii#) 



gi{R, Zk) + Rd R gi(R, z k ) 



1 

2na J 
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Va 2 - R 2 Jo 
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2na 



g 2 (R, Zk) - a 
RdR 
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1 
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y/a 2 - R 2 Jo 



d9 cos(fi2&) 
g 3 (R, z k ) + Rd R g 2 {R, z k ) 

2tt 



X 



1 

2na 



ga(R,Zk) - a 
RdR 
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1 



M2 

R 



d6 cos(;iii 



g 2 (R, z k ) + Rd R g 3 (R, z k ) 



2n 



X 



9i{R, Zk) - a 



Va 2 - R 2 Jo 
Mi" 



d R - 



R 



d6 cos(fi 2 9) 
gi(R, z k ) + Rd R gi{R, z k ) 



(C.l) 



(C.2) 



(C.3) 



(C.4) 



When /Xi,/X2 7^ 0, then tp^>(t n+ i, 0, z k ) — by virtue of the angular integration, in agreement with the 
boundary conditions in Eqs. ( |l~9| ) to (22). When j z = | (the case j z — —\ is similar), the preceding 
equations become: 



i){j(t n+ i,0,Zk) 



RdR—= 

o Va 2 - R? 



31,3 {R, z k ) - a 



Or 



g2,i(R, z k ) + Rd R g it3 (R, z k ) 



(C.5) 



while tp^li^n+i, 0, Zk) = 0. Using the same strategy as above, the cubic Hermite interpolation becomes 

g h ',i(r,z k ) = g h ,AU,k)S 1 (ar))+g h 'A^,k)S 2 (ttr)) + r{ h ' ) m h , A (l,k)S 4 (ar)), (C.6) 

g h 'Ar,z k ) = fr,2(U)S 2 KH) + rf ) m k ', 2 (0,i)S 3 KW) + r f ) m k - 2 (U)S 4 K(r)), (C.7) 

9h',a(r,z k ) = ffv,3(0,fc)Si«(r))+flv,3(l,*:)S , 2 (f(r)) + r[ V) mv,a(l,fc)S 4 (€(r)), (C.8) 

gv A (r,z k ) = ^,4(l,fc)^(^r)) + r^ ) m^, 4 (0,fc)53(^(r)) + r^' ) m^,4(l,A ; )54(e(r)), (C.9) 
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where £(r) := and where we set m^i 13(0, k) = according to the boundary conditions in Eqs. 
22 Also, the remaining derivatives in r — are approximated by 
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to 



^',2,4(0, k) 



9h'(l,k) 



2r 



(/»' 



(CIO) 



where we used the fact that the wave function components 1 and 4 are antisymmetric around r — 0. Thus, 
with these techniques, the boundary conditions at r — is fulfilled exactly by the interpolating polynomial. 



The next step is the substitution of the interpolant into Eqs. (C.l) to (C.4), resulting in integrals of the 
form 



dR , = 
s/a 2 - R 2 



:V7T r (l + k) 



(C.ll) 



2 r(i + i)- 

It is then possible to evolve the wave function at r — by using points within the domain, while preserving 
the boundary conditions. Also, it is clear that we circumvent the coordinate singularity problem as the 



integration in Eq. (C.ll) is well-defined. 



Appendix D. Solution of 2-D wave packet 

In this Appendix, the analytical solution for the time evolution of a 2-D free wave packet with azimuthal 
symmetry is computed. The solution can be computed in polar coordinates by considering the solution in 
Cartesian coordinates. We consider an initial wave packet for a massive spin-up electron at rest. The initial 
wave function in Cartesian coordinates is given by 



Ht = 0,x,y) =N 



e i(^) 2 



and its Fourier transform by 



ip(t = 0,p x ,p v ) = 4ttA 2 U 



The 2-D Dirac equation we want to solve is given by 

id t Tp(t,p x ,p y ) = [ca x p x + ca y p y + /3mc 2 ] ip(t,p x ,p y ), 
here expressed in Fourier space. The solution to this equation is then simply 

= I 4 cos (Et) - i ? ; ' : — .m! El . 

xtKO,^,^) 



ca x Px + cotyPy + fimc 2 
E 



(D.l) 



(D.2) 



(D.3) 



(D.4) 



(D.5) 
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where E = yp 2 c 2 + p 2 c? + m 2 c A . This last equation can be Fourier transformed back to real space. Then, 
using polar coordinates (p 2 = p 2 +Py) and properties of Bessel functions, we get the solution as 

poo 

Mt,r,z) = 27VA 2 / dppe- A2p2 J (pr) 



ii 



x 



. TOC 2 



cos(Et) — i sin(Et) 

E 



(D.6) 

Mt,r,z) = (D.7) 

Mt,r,z) = (D.8) 

Mt,r,z) = 2JVA 2 e lS / dppJ^pr)^ sm(M) e - AV (D.9) 

Jo E 

where J n {z) is the Bessel function of the first kind. Looking at the angular dependence, we see that this 
solution corresponds to a wave function with j z = 1/2. 

There is no known analytical form for these integrals in the general case [55] . However, using high 
accuracy numerical integration, this result can be used to validate our numerical method and to analyze the 
operator splitting. 
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